Liquid-liquid phase transitions in supercooled water studied by computer simulations 

of various water models 
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Liquid-liquid and liquid-vapor coexistence regions of various water models were determined by MC 
simulations of isotherms of density fluctuation restricted systems and by Gibbs ensemble MC sim- 
ulations. All studied water models show multiple liquid-liquid phase transitions in the supercooled 
region: we observe two transitions of the TIP4P, TIP5P and SPCE model and three transitions of 
the ST2 model. The location of these phase transitions with respect to the liquid-vapor coexistence 
curve and the glass temperature is highly sensitive to the water model and its implementation. We 
suggest, that the apparent thermodynamic singularity of real liquid water in the supercooled region 
at about 228 K is caused by an approach to the spinodal of the first (lowest density) liquid-liquid 
phase transition. The well known density maximum of liquid water at 277 K is related to the second 
liquid-liquid phase transition, which is located at positive pressures with a critical point close to the 
maximum. A possible order parameter and the universality class of liquid-liquid phase transitions 
in one-component fluids is discussed. 
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I. INTRODUCTION 

Some one-component substances show several amor- 
phous states at low temperatures. This sug- 
gests the possible existence of phase transitions be- 
tween them. Evidence for such phase transi- 
tions in supercooled water was obtained both from 
experiments^* 2 *^ and simulations£*§iL2i2ii2iiiii2 (see also 
recent reviewsi^i^i^i^) . The presence of phase tran- 
sitions between different liquid (glassy) states of an 
isotropic fluid in the supercooled state (below the freez- 
ing temperature) may noticeably influence its properties 
above the freezing temperature. In particular, at low 
temperatures an anomalous, singularity-like behavior of 
some thermodynamic properties is observed in several 
fluids^ including water pi&iiLSI Such behavior can be ex- 
plained by the influence of a distant liquid-liquid phase 
transition with a critical point located in the supercooled 
region. 5 

Liquid-liquid phase transitions of a one-component 
fluid may occur in a thermodynamic region where one or 
both liquid phases are metastable with respect to crystal- 
lization or evaporation. This makes experimental studies 
of such phase transitions difficult or even impossible. For- 
tunately, this is not the case for computer simulations, 
which can be used to study fluids in metastable or even 
unstable states. In principle, coexistence of two liquid 
phases with an explicit interface between them can be 
obtained by direct simulations in the constant-volume 
NVT ensemble. However, this demands the use of large 
systems, which can not be equilibrated at low tempera- 
tures during reasonable computer simulation times. De- 
creasing the system size causes a narrowing of the den- 
sity range, where the coexistence of two phases with ex- 
plicit interface can be obtainedi^i In a small enough sys- 
tem, whose size does not exceed some "threshold" value, 
a phase separation with a persistent interface can be 
prevented completely and a continuous isotherm, join- 
ing the two states at subcritical temperatures, can be 



simulated. The presence of a van der Waals-like loop di- 
rectly evidences the sub-critical character of an isotherm. 
This allowed to detect a liquid-liquid phase transition in 
ST2 5 i 6 i 9 i 22 i 23 i 24 and TIP5B& water models. Note, that 
the unavoidable tendency of the simulated metastable 
and unstable systems to phase separation ultimately dis- 
torts the computed isotherms; 25 - Moreover, the threshold 
value of the system size, which prevents phase separa- 
tion in the simulation box in the whole density range of 
the two-phase coexistence region, can not be estimated 
in advance. These factors represent intrinsic limitations 
of the applicability of simulations in the NVT ensemble 
to study phase transitions, which can not be overcome by 
varying the system size or by improving the sampling. In 
particular, these shortcomings can lead to the absence of 
a van der Waals-like loop in a sub-critical isotherm close 
to the critical temperature, resulting in a lower value of 
the estimated critical temperature of the phase transi- 
tion. 

Simulations in the constant-pressure NPT ensemble 
can also be used to study low-temperature liquid-liquid 
phase transitions. An abrupt change of the density along 
an isobar at some temperature indicates a phase tran- 
sition. However, it is difficult to distinguish a phase 
transition from a sharp but continuous change of the 
densityi2£i2£ Only stable and metastable but not me- 
chanically unstable states can be explored in the NPT 
ensemble and, therefore, a continuous isotherm, joining 
two phase states at subcritical temperatures, can not 
be obtained. Simulations of subcritical isotherms in the 
NPT ensemble are accompanied by hysteresis phenom- 
ena, which may serve as an indicator of a first-order 
phase transition. However, the extension of a hystere- 
sis loop depends on the system size and simulation time. 
As in the case of the NVT ensemble, the extension of 
such simulations beyond the stable states could hardly 
be controlled. 

The study of phase transitions by simulations in canon- 
ical ensembles can be substantially improved by forcing 
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the system to remain homogeneous even in metastablc 
and unstable states. This can be achieved by restricting 
the density fluctuations on a mesoscopic scalei 2 ^ 2 ? Sim- 
ulations of isotherms in such restricted ensembles enable 
the location of liquid-liquid phase transitions at super- 
cooled temperatures4fl*Ai In the present paper we use this 
method to find liquid-liquid phase transitions of various 
water models in the supercooled region. In some cases the 
densities of the coexisting liquid phases were determined 
directly by simulations in the Gibbs ensemble To 
locate liquid-liquid transitions with respect to the liquid- 
vapor coexistence curve, the latter was simulated in the 
Gibbs ensemble down to low temperatures. The obtained 
phase diagrams are used to analyze the sensitivity of the 
liquid-liquid and liquid- vapor phase transition to the pe- 
culiarities of a given water model and its implementation. 
From these data the possible location of the liquid-liquid 
phase transitions in real water is inferred. 



II. SIMULATION METHODS 

ST2,22 TIP4P,2a TIP5R2I and SPCE 32 water models 
were used in our studies. In all cases a simple spheri- 
cal cutoff R c of the intermolecular interactions was used 
with long-range corrections for the Lennard- Jones po- 
tentials (LRLJ corrections). The value of R c was 9 A 
in all cases, except the TIP4P model, where R c was 
8.5 A in the simulations of the liquid-vapor coexistence 
curve. In accord with the original parametrization of 
these models^2iiffii2ii2£ no long-range corrections for the 
Coulombic interaction (LRCI corrections) were included. 
Additionally, we used the ST2 water model including a 
reaction field to account for LRCI corrections (ST2RF 
water model hereafter), which has already been stud- 
ied intensively in the supercooled regioni^S*2i2SiSiS4 The 
liquid-vapor coexistence curve of various water models 
was obtained by Monte Carlo (MC) simulations in the 
Gibbs ensemble (GEMC)^ The coexistence of different 
water phases at supercooled temperatures was studied 
by GEMC simulations and by simulations of isotherms 
in the density fluctuation restricted APT ensemblei 25 i 28 



A. Gibbs ensemble MC simulations of the 
liquid-vapor and liquid-liquid coexistence 

The liquid-vapor coexistence of water was simulated 
by direct equilibration of the two coexisting phases in 
the Gibbs ensemble^ The total number A of molecules 
in the two simulation boxes was about 600, except for 
the case of the TIP4P model, where A was about 300 
to 400, depending on temperature. The use of efficient 
techniques for the molecular transfers (more details can 
be found elsewhere^) allowed to extend the simulated 
coexistence curves deeply into the supercooled region. 
Even at extremely low temperatures the number of suc- 
cessful molecular transfers in the course of the simula- 



tion runs always exceeded 10A and achieved 100A at 
high temperatures. The obtained liquid densities did not 
depend on the initial configuration of the simulated sys- 
tem even at the lowest studied temperatures. This was 
achieved by efficient sampling of the system due to molec- 
ular transfers and also by the absence of phase separation 
in the simulation box, which is intrinsically avoided in the 
Gibbs ensemble. The critical temperature and critical 
density of the liquid-vapor coexistence curves were esti- 
mated by fitting the order parameter with a scaling law 
including one correction term and by fitting the diameter 
to a second order polynomial. The coexisting densities 
above the temperature of the liquid density maximum 
were used for this fitting procedure. 

The coexistence between liquid water phases of differ- 
ent densities was also studied by GEMC simulations. The 
extremely low acceptance probability for molecular trans- 
fers between two dense liquid phases practically prevents 
such simulations at temperatures below ^260 K and den- 
sities above ~1.3 g cm -3 . Another complication is caused 
by the possible small difference between the densities of 
the coexisting liquid phases. GEMC simulations develop 
two coexisting phases, when the average density p av of 
the two simulation boxes is in the two-phase region. The 
initial choice of this average density determines the box 
sizes and numbers of molecules in both boxes, which have 
to be adequate to reproduce the properties of both coex- 
isting phases. This essentially restricts the range of p av 
, which are appropriate for such studies. Therefore, we 
have equilibrated several initial configurations of liquid 
water with different densities in the interval from 0.90 to 
1.25 g cm" 3 by MC simulations in the NVT ensemble. 
Then, several pairs of these configurations with density 
differences from 0.10 to 0.20 g cm -3 were used as start- 
ing configurations for GEMC simulations. The number 
of successful molecular transfers was about A at T = 
260 K, about 2A at T = 270 K and essentially more at 
higher temperatures. Any choice of the average density 
p av at T > 270 K resulted in comparatively fast (less one 
month of computer time on a GHz processor) equaliza- 
tion of the densities in the two boxes. The same result 
was obtained also at T = 260 or 270 K, when p av > 
1.10 g cm" 3 . However, when the average density p av was 
about ^0.95 and ~1.05 g cm" 3 the water densities in the 
two boxes did not tend to equalize at T = 260 and 270 
K even during extremely long simulation runs (several 
months of computer time on a GHz processor) but a sta- 
ble liquid-liquid equilibrium was reached. It should also 
be noticed, that we have failed to get the coexistence of 
two liquid water phases, when starting from two equal 
densities in the two simulation boxes. 



B. MC simulations in the restricted NPT ensemble 

An extension of the isotherm into metastable or un- 
stable regions can be achieved by simulations in the den- 
sity fluctuation restricted ensemble which allows the 
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FIG. 1: Liquid-vapor coexistence curve of the ST2 and 
ST2RF water models (circles). Experimental data&2& are 
shown by the solid lines. 

simulation of an artificial single-phase state. The simula- 
tion box with N molecules has to be divided into n equal 
subcells. In the general case, the number N of molecules 
in the i-th subcell should satisfy the following constraint: 



N/n -5N<N,< N/n + SN, 



(1) 



where SN is the maximum deviation of the number of 
molecules in the subcell from the average value N/n, al- 
lowed in simulations. Such a constraint can be fulfilled in 
MC simulations by rejecting those moves, which violate 
Eqffl The value N/n should be small enough to exclude 
phase separation in each subcell, whereas too small val- 
ues of N/n could result in the freezing of the simulated 
liquid. MC simulations in the restricted NVT ensemble 
were successfully applied for reproducing the liquid- vapor 
coexistence curve of a LJ fluidiS 5 ^ 2 ^ In the case of liquid- 
liquid coexistence in a one-component fluid the densities 
of the coexisting phases could have rather close values 
(difference < 10 %). So, even SN = 1 could result in a 
phase separation in a simulation box with several hun- 
dreds molecules. Therefore, we imposed SN = 0, which 
could noticeably suppress the density fluctuations in the 
NVT ensemble. 25 - 28 Imposing this local restriction in an 
NPT simulation, the fluctuation of the average density 
is not restricted. This method can be applied for sim- 
ulations of stable and metastable, but not for unstable 
states. 

In our implementation of the fluctuation restricted 
NPT ensemble, the cubic simulation box with 513 water 
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FIG. 2: Liquid densities of ST2 (open circles) and ST2RF 
(solid diamonds) water models along the liquid-vapor coexis- 
tence curve. Simulation data at P ~ for the ST2 model 29 
and the ST2RF model- are shown by solid circles and stars, 
respectively. Experimental data;— shifted upwards by 28° , 
are shown by the solid line. 



molecules was divided into 27 equally sized cubic sub- 
cells which contain an equal number (19) of molecules. 
The number of molecules in each subcell was kept un- 
changed in the course of the MC simulations. Typically 
up to 2*10 5 molecular moves per molecule were done in 
the course of a MC simulation run. Several isotherms 
between 150 and 300 K were simulated for each studied 
water model. To test the effect of the restrictions on the 
thermodynamic properties of the model water, we simu- 
lated the isobar P — in the restricted NPT ensemble 
for the ST2 model from 200 to 450 K. The liquid den- 
sities obtained by this method were compared with the 
saturated liquid densities, obtained by the GEMC simu- 
lations (see Fig. 2). Additionally, the isotherm T — 260 K 
of the ST2 model was simulated in a simulation box with 
621 molecules, divided into equally sized cubic subcells 
which contain 23 molecules each. 

Since neither the location of the liquid-liquid phase 
transitions (except for one transition of the ST2RF 
mo ^ e i5 i 6 4 9 4 22 i 23 4 24^ nor even their number are known a 

priori, we explored in detail a wide interval of density 
(from ^0.8 to ~1.4 g cm~ 3 ) and pressure (from ~-6 
to ^+10 kbar). To detect all branches of the studied 
isotherm, we performed simulations in the following way. 
Several configurations of liquid water with different den- 
sities were equilibrated in the restricted NVT ensemble. 
Then, each prepared configuration was used as an initial 
one for simulations in the restricted NPT ensemble at 
several chosen pressures. The systems converge to one 
(at supercritical temperatures) or more (at subcritical 
temperatures) densities, which correspond to all possi- 
ble stable and metastable states of liquid water at the 
considered pressure and temperature. At low tempera- 
tures such convergence demanded extremely long Simula- 
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FIG. 3: Liquid- vapor coexistence curve of TIP4P and TIP5P 
water models (circles) . Simulation data of Refl37l and Refl3H 
are shown by triangles and stars, respectively. Experimental 
data&— are shown by solid lines. 



tion runs (several months in some cases). Subsequently, 
the density along the isotherms in each detected water 
phase was studied by increasing (decreasing) the pres- 
sure in small steps (0.1 kbar in some cases). At some 
pressure, which roughly corresponds to the stability limit 
(spinodal), the density of a considered water phase in 
the course of a simulation run rapidly shifts to another 
branch of the isotherm, corresponding to another wa- 
ter phase. This approach should allow for discovering 
all branches of a water isotherm, whose densities at the 
same pressure differ more than about 0.03 g cm~ 3 . This 
value approximately corresponds to the thermodynami- 
cal fluctuations of the density in our finite system. Note, 
that we checked the influence of the initial conditions of 
the simulation runs. Even at the lowest temperatures 
we observed convergence to the given results. This was 
due to the extremely long simulation runs and the ab- 
sence of phase separation in the simulation box, which is 
intrinsically avoided in the restricted ensemble. 

The obtained isotherms do not allow to locate accu- 
rately the coexistence pressure and the densities of the 
coexisting phases due to the missing parts in the region 
of the unstable states. Moreover, the density is not nec- 
essarily the appropriate order parameter to describe the 
liquid-liquid phase transition of a one-component system. 
Therefore, we attributed the coexistence pressure to the 
average of the maximal and minimal pressures of the 



TABLE I: The values of the liquid- vapor critical temperature 
T c and critical density p c , temperature of the liquid density 
maximum T max at the liquid-vapor coexistence curve or at 
zero pressure, temperature of the glass transition T g , obtained 
from simulations of the studied water models. Experimental 
data are also shown. 



Model 


T c 


Pc 


Tmax 


Tc T m a x 


T g 
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(K) 


ST2 


550.2 


0.286 


305 


246 


235 








~30C 29 






ST2RF 


536.2 


0.286 


320 


216 


255 








~33022 






TIP4P 


580 


0.328 


~ 250 


330 


180 




57S 38 




263^ 








588. 2 41 


0.315^i 








TIP5P 


537.7 


0.290 


~ 275 


263 


215 




521.3 41 


0.337^ 


275 31 








546 38 










SPCE 


623.3 


0.319 






220 




64C 42 


0.29^ 


235^ 




17?S 




630^ 


0.295^ 


255i& 




~21C 45 




638.C 46 


0.273^ 


~24C 39 




188 47 


Exp. 


64? 3 £ 


0.326^ 


277.1^ 


374.6 


~16C 48 












~136 15 



overlap region, where the two branches have two distinct 
densities at the same pressure. Above the glass transi- 
tion temperature T g , the accuracy of such estimations is 
comparable with the accuracy of the GEMC simulation 
data. Far below T g , the density ranges of the different 
branches of the isotherm overlap so strongly, that the 
densities of the coexisting phases can be estimated only 
very roughly. Note, that at T < T g , we do not infer the 
existence of some phase transition from abrupt variations 
of the density in a small pressure interval. Additional in- 
formation about the sign of the coexistence pressure we 
can obtain by attributing the density of the liquid branch 
of the liquid- vapor coexistence curve at the same temper- 
ature to zero pressure. 



III. RESULTS 

A. Liquid-vapor coexistence of water from the 
critical point to the glassy state 

1. Liquid density along the liquid-vapor coexistence curve 

The liquid-vapor coexistence curves of the ST2 and 
ST2RF water models are shown in Fig. 1 . At high temper- 
ature in both models the liquid densities are essentially 
lower than the experimental values due to the strongly 
underestimated critical temperatures T c (see Table I). 
At ambient temperatures the liquid density of the ST2 
model is close to the density of real water, whereas in 
the ST2RF model it remains significantly underestimated 
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FIG. 4: Liquid densities of TIP4P and TIP5P water mod- 
els along the liquid- vapor coexistence curve (open circles). 
Simulation data at P ~ for the TIP4P model: squares, 00 
triangles, 6 stars, 51 diamonds^- Simulation data at P ~ 
for the TIP5P model: stars;— diamondsi— Experimental 
data— ^ are shown by solid lines. 
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FIG. 5: Liquid-vapor coexistence curve of SPCE water model 
(open circles). Simulation data from Refli^. Refl44l and 
R.eflH^ are shown by dashed line, solid triangles and solid 
squares, respectively. Experimental data—*— are shown by 
solid line. 



(Fig.2). Our results of the GEMC simulations (Fig.2) 
agree well with the results of molecular dynamics (MD) 
simulations of the ST22S and ST2RB& models at P w 0. 
The ST2 model shows the density maximum at about 305 
K and well reproduces the experimental liquid density in 
a temperature range of about 70°, when the latter is 
shifted upwards by 28° (see Fig.2). At all temperatures 
the liquid density of the ST2RF model is significantly 
below the experimental values. In the ST2RF model the 
temperature T max of the density maximum (about 320 
K) is slightly higher than in the ST2 water, and the ab- 
solute value of the liquid density pi is about 7% below 
the experimental value. The experimental density pi of 
about 1 g cm~ 3 can be reproduced in the ST2RF model 
by applying a positive pressure of about 800 bars^Si 
The drastic density difference between ST2 and ST2RF 
models is caused by the application of the reaction field 
method for the treatment of the LRCI corrections in the 
latter case, which apparently strengthens the water hy- 
drogen bonding. The shift of the liquid density in the 
ST2RF model agrees with the predictions of the modi- 
fied van der Waals model for water — 1 The characteristic 
temperature interval from T c to T max is only 66% of ex- 
perimental value for the ST2 model and only 58% for 
the ST2RF model (Table I). Note, that with decreasing 
temperature a density minimum follows the density max- 
imum in both models (see Figs.l and 2). 

The striking feature of the liquid-vapor coexistence 
curve of the ST2 model is a break of the liquid branch at 
T = 270 K (Figs.l and 2). At this temperature two liq- 
uid phases with different densities (about 0.91 and 0.97 g 
cm -3 ) coexist with the vapor phase. This is not the case 
for any other studied temperatures, including the nearest 
ones T = 265 K and T = 275 K, where only one liquid 
phase coexists with the vapor. Evidently, there is a triple 
point at T = 270 K, where two first order phase tran- 



sitions (liquid- vapor and liquid-liquid) meet each other. 
The gradual change of pi in the whole studies temper- 
ature range (Fig.2) indicates the absence of the triple 
point in the ST2RF model. The disappearance of the 
triple point with the strengthening of the water hydro- 
gen bonding qualitatively agrees with the predictions of 
the modified van der Waals model for water — i 

The liquid- vapor coexistence curves of the TIP4P and 
TIP5P water models are compared with the experimen- 
tal data in Fig. 3. Upon cooling the liquid density in 
both models passes through a maximum and afterward 
through a minimum. Both models underestimate the 
critical temperature (Table I), that essentially depresses 
the liquid density at high temperatures, especially in the 
case of the TIP5P model. For the TIP4P water model 
the characteristic temperature interval from T c to T max 
is about 88% of the experimental value, whereas for the 
TIP5P model it is only 70%, which is close to the value 
for another five-site water model, the ST2 model (see 
Table I). Most of the data points for the TIP4P model 
are from our previous paper^ where comparisons with 
other simulation studies of the coexistence curve can be 
found (see Fig. 5 in Ref l34T) . Depending on the imple- 
mentation of the TIP4P and TIP5P water models (long- 
range corrections for intermolecular interactions, cutoff, 
system size, etc.), the density of the saturated liquid at 
room temperature varies within 1-2 %. In particular, us- 
ing the LRCI corrections causes a slight decrease of the 
liquid density, whereas including the LRLJ corrections 
causes the opposite trend. 

The details of the particular implementation of the wa- 
ter model and the employed simulation method gain more 
importance with decreasing temperature (Fig. 4). At am- 
bient temperatures the liquid density of the TIP4P wa- 
ter model is not very sensitive to the simulation details. 
However, two qualitatively different courses of the liquid 
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FIG. 6: Liquid densities of the SPCE water models along 
the liquid- vapor coexistence curve (open circles). Simulation 
data at P ~ 0: triangles up^~ triangles down;— squares;— 
stars— and diamonds. 56 Experimental data*-*— are shown by 
the solid line. 



density are obtained in the supercooled region: in our 
simulation pi starts to increase with decreasing temper- 
ature at T < 225 K, whereas in other simulations^^ it 
continues to decrease. Apparently, this difference origi- 
nate from the use of LRLJ corrections in the former case 
and LRCI corrections in the latter case. 

The TIP5P water model seems to be even more sensi- 
tive to the simulation details In the region of the 
density maximum the use of LRCI corrections 3 ^ causes 
a decrease of pi by about 2%, whereas the use of LRLJ 
corrections (our simulations) causes a slight increase of 
pi in comparison with the original TIP5P model, 31 which 
was parameterized without any long-range corrections for 
intermolecular interactions (Fig. 4). In the supercooled 
region the increase of pi due to the LRLJ corrections 
achieves about 3%. 

Among the considered water models, the SPCE model 
gives the most realistic value of T c (see Table I), and ac- 
cordingly the closeness of the simulated and experimen- 
tal values of the coexisting densities at high temperatures 
(Fig. 5). The characteristic temperature interval from T c 
to T max is also close to the experimental value (Table 
I). Application of the LRCI corrections2244i£2 depresses 
the liquid density by about 2-3% similarly to other water 
models. At low temperatures the behavior of the den- 
sity of the saturated liquid in the SPCE model seems to 
be more sensitive to the simulation details in comparison 
with other models (Fig. 6). In the region of the density 
maximum the value of pi is the largest, when only LRLJ 
corrections are taken into account (Fig. 6, circles). It sub- 
sequently decreases, when such corrections are excluded 
(Fig. 6, solid triangles), when both kinds of the long-range 
corrections are included (Fig. 6, diamonds) and when only 
LRCI corrections are included (Fig. 6, squares). The use 
of a smaller spherical cutoff results in the lowest liquid 
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FIG. 7: Potential energy U in the liquid phase along the 
liquid-vapor coexistence curve of the ST2 and ST2RF water 
models (open circles) . Data points in the vicinity of the liquid- 
liquid transition of the ST2 model at T ~ 270 K are shown 
in the inset. The values of U, obtained from simulations of 
the ST2RF model at P w Of 6 - are shown by stars. 



density (Fig. 6, open triangles). 

Some features of the temperature dependence of the 
saturated liquid density at ambient and supercooled tem- 
peratures seem to be general for all studied water models, 
despite the high sensitivity of the phase diagram to the 
model and its implementation. All studied water mod- 
els show the well-known density maximum (see Table I) . 
This maximum is followed by a density minimum with 
decreasing temperature. Previously, we observed the liq- 
uid density minimum along the liquid-vapor coexistence 
curve for TIP4P water 34 i 57 i 58 (Figs.3 and 4), ST2 waterifi 
(Figs.l and 2) and TIP5P watenii (Figs.3 and 4). Now, 
the ST2RF water model also shows a pronounced den- 
sity minimum (Figs. 1 and 2). This well agrees with the 
recent observation of a density minimum of this model at 
negative pressures— ^ Although we are not able to locate 
both a density minimum and a maximum in our simula- 
tions of SPCE water (Fig. 5), the existence of the density 
minimum seems to be quite possible for other implemen- 
tations of this model (Fig. 6). Note, that the liquid den- 
sity minimum around 185-195 K was also observed upon 
temperature quench of a water model, proposed by Guil- 
lot and Guissani— i 
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FIG. 8: Potential energy U in the liquid phase along the 
liquid- vapor coexistence curve of the TIP4P and TIP5P water 
models (open circles) . The values of U, obtained from simula- 
tions at P ~ are shown by squares;— stars- and trianglesi— 



2. Potential energy of liquid water along the liquid-vapor 
coexistence curve and glass transition 

The temperature dependence of the potential energy 
U of the water molecules talong the liquid-vapor coex- 
istence curve is shown in Fig. 7 for the ST2 model. The 
discontinuity of U (of about 0.23 kcal/mol) at T = 270 
K (see inset in Fig. 7) reflects the liquid- liquid phase tran- 
sition at the triple point. The slope of U(T) does not 
change noticeably when crossing the triple point, indicat- 
ing the essentially liquid nature of the coexisting phases. 
A slower change of U(T), typical for the glassy state, is 
observed at low temperatures. This allows an estimate of 
the glass transition temperature of the ST2 water model 
at zero pressure as T g w 235 K. 

The temperature dependence of the potential energy 
U of the STRF water model indicates a transition to the 
glassy state at T g w 255 K, i.e. slightly higher then in 
the ST2 model. Two different slopes of U (T) above T g 
evidences two quite different kinds of liquid water. The 
continuous transition between these two liquids occurs at 
about 300 K, i.e. about 20° below the temperature of the 
density maximum (Fig. 7, Table I). 

Transitions to glassy states are also seen for the TIP4P 
and TIP5P water models (Fig. 8). The estimated values 
of the glass transition temperature are: T g w 180 K for 
the TIP4P model and T g w 215 K for the TIP5P model. 
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FIG. 9: Potential energy U in the liquid phase along the 
liquid-vapor coexistence curve of the SPCE water model 
(open circles). Values of U, obtained from simulations at P 
ps 0: triangles up;— triangles down~ squares;— stars— and 
diamonds^ 

Note, that the use of the LRCI correction (Fig. 8, squares 
and stars in the upper panel and stars in the lower panel) 
results in noticeably lower potential energies at T < 225 
K for the TIP4P model and T < 275 K for the TIP5P 
model. This is obviously connected with the essentially 
different densities of water in the two implementations of 
the TIP4P and TIP5P water models in these tempera- 
ture ranges (see Fig. 4) . The glass transition temperature 
of the SPCE water model noticeably depends on its im- 
plementation (see Table I) and from our simulation study 
we estimate glass transition temperature T g « 220 K for 
this model (Fig. 9). 

3. Heat capacity along the liquid-vapor coexistence curve 

The heat capacity at constant pressure C p is the tem- 
perature derivative of the enthalpy and relates to the 
potential energy of a system through the equation: 

c,= ( a ^ P 9 v T +ZRT) ) p < 2 > 

Along the liquid-vapor coexistence curve, the contribu- 
tion of P{d V/9 T)p to C p is negligible up to tempera- 
tures of about 450 to 500 K. We have computed C p from 
2- or 3-point running averages of the enthalpy U + 3RT 
followed by central differentiation. The data for the low- 
est temperature was estimated using the forward formula 
for the derivative. 

Our simulation results are compared with the experi- 
mental data and other simulation studies of C p in Figs. 10 
and 11. The accuracy of the obtained C p data directly 
depends on the number of simulated temperature points. 
The available simulation data allow an estimation of the 
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FIG. 10: Heat capacity of liquid water along the liquid- 
vapor coexistence curve: our simulation (open circles) and 
experiment—*— (thick lines). Experimental data are shifted 
by +35° and +55° for comparison with the ST2 and ST2RF 
water models, respectively. Dotted line shows an estimate of 
C p , when the discontinuity at the liquid-liquid transition is 
not taken into account. Dashed lines indicate the values of 
C p , estimated from the linear approximations of U, shown by 
solid lines in Fig. 7. 



FIG. 11: Heat capacity of liquid water along the liquid- 
vapor coexistence curve: our simulation (open circles) and 
experiment— *— (thick lines). Data from other simulation 
studies at P ~ 0: TIP4P model— (squares, upper panel), 
TIP4P-Ew model 54 (stars, upper panel), TIP5P model 31 
(squares, middle panel), TIP5P-Ew model— (stars, middle 
panel). Dashed lines indicate the values of C p , estimated 
from the linear approximations of U, shown by solid lines in 
Figs.8 and 9. 



average value of C p over a wide temperature interval, 
whereas the temperature dependence of C p can be ana- 
lyzed only qualitatively. In some temperature intervals 
the temperature dependence of U is also approximated 
by a linear behavior (see solid lines in Figs. 7 to 9) and 
the corresponding average values of C p are shown by hor- 
izontal dashed lines in Figs. 10 and 11. The value of C p 
in the high-temperature liquid is about twice the C p in 
the glassy state. 

Different temperature dependences of C p are observed 
for the different water models. An increase of the heat 
capacity of the ST2 water model is clearly seen with de- 
creasing temperature down to 270 K, i.e. to the triple 
point of the liquid-vapor and liquid-liquid phase transi- 
tion (Fig. 10, upper panel). An estimation of C p , which 
neglects the discontinuous drop of the potential energy 
U at T\ = 270 K (see inset of Fig. 7), results in a much 
stronger apparent increase of the heat capacity (dotted 
line in Fig. 10). Such an artificially strong increase of the 
estimated value of C p could appear, if the liquid-liquid 
phase transition is smeared our due to the limitations 
of the simulation method. The heat capacity of the low 
temperature liquid phase can hardly be estimated due to 
the small temperature interval between the liquid-liquid 
phase transition and the glass transition (see Fig. 7), but 
seems to be of the same magnitude as above the triple 
point. 

The ST2RF model also shows an increase of C p upon 
cooling (Fig. 10, lower panel). At about 300 K, where a 
gradual transition to the low density form of liquid water 



occurs, C p increases abruptly to a value of about 0.04 
kcal mol -1 K _1 , and remains at such high values till the 
transition to the glassy state. So, the heat capacity of 
the second (low-density, low-temperature) liquid in the 
ST2RF model is almost two times higher, than the heat 
capacity of the ordinary liquid, which coexists with the 
vapor in a wide temperature range up to the liquid- vapor 
critical point. 

Similarly to the ST2 and ST2RF models, the TIP5P 
water model shows a noticeable increase of C p with de- 
creasing temperature (Fig.ll, middle panel). This be- 
havior is interrupted by the transition to a glassy state. 
The absolute values of C p in the TIP4P and SPCE wa- 
ter models are rather close to the experimental data in a 
wide temperature range (Fig.ll). However, these models 
do not show the noticeable increase of C p upon cooling, 
observed in real water. 



B. Liquid-liquid phase transitions of water at low 
temperatures 

1. ST2 model 

Isotherms of the ST2 water model, obtained by MC 
simulations in the restricted NPT ensemble at T — 235 
to 290 K, are shown in Fig. 12. No liquid-liquid transi- 
tions could be detected at T = 290 K, where the liquid 
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density gradually changes with pressure. At T = 275 
K a liquid-liquid coexistence appears at slightly negative 
pressure (P w -0.2 kbar), the densities of the coexisting 
liquid phases are about ~0.94 and ~0.98 g cm -3 . At 
T = 260 K this transition shifts to positive pressure (P 
ss +0.4 kbar) and the two-phase region becomes wider 
(from about 0.90 to 1.01 g cm -3 ). So, at some temper- 
ature between 275 K and 260 K this liquid-liquid transi- 
tion should pass through zero pressure, i.e. should meet 
the liquid- vapor phase transition at a triple point. In- 
deed, such a triple point at T = 270 K was found at the 
liquid-vapor coexistence curve of the ST2 model (Figs.l 
and 2). The densities of the coexisting liquid phases at 
the triple point, estimated from the GEMC simulations 
of the liquid-vapor coexistence, ~0.91 and ~0.97 g cm -3 , 
are in good agreement with the MC simulations in the re- 
stricted APTensemble. Moreover, extensive direct equi- 
libration of the two liquid phases in the Gibbs ensemble 
also gives a rather similar coexistence region: between 
0.92 and 1.00 g cm" 3 at T = 270 K and between 0.90 
and 1.02 g cm" 3 at T = 260 K (see Fig.4 in RefEJ. 
Pressure and densities of these coexisting phases do not 
change noticeably upon cooling to T = 235 K. 

At T — 260 K the ST2 water model shows a second 
liquid-liquid phase transition, occurring at positive pres- 
sure (P « +0.5 kbar) in the density interval from ~1.02 
to ~1.11 g cm~ 3 . This interval agrees well with the densi- 
ties of the coexisting phases, obtained directly in GEMC 
simulations: 1.03 and 1.12 g cm" 3 at T = 270 K and 
1.04 and 1.12 g cm" 3 at T = 260 K (Fig.4 in RefEjJ. 
The densities of the coexisting phases do not change no- 
ticeably at T — 235 K, whereas the coexistence pressure 
slightly increases up to P « +0.85 kbar. A third liquid- 
liquid phase transition of the ST2 water model appears 
at T = 235 K. It is located at positive pressure (P w 
+0.95 kbar) (Fig. 12) and the densities of the coexisting 
liquid phases are ~1.12 and ~1.19 g cm -3 . 

The location of the thus observed liquid-liquid coex- 
istence regions with respect to the liquid branch of the 
liquid-vapor coexistence curve is shown schematically in 
upper panel of Fig. 13 (see also Fig.4 of RefJiOi. The 
simulation results evidence, that the first (the lowest- 
density) liquid-liquid phase transition of the ST2 model 
should have a critical point of demixing at slightly neg- 
ative pressures with a critical temperature between 275 
K and 290 K and a critical density between ~0.94 and 
^0.98 g cm~ 3 . It seems likely, that the second and 
third phase transitions also end at corresponding criti- 
cal points. However, the observed temperature evolution 
of the coexisting densities of these transitions does not al- 
low to exclude the possibility of a common critical point. 
If so, the narrow one-phase region near 1.1 g cm~ 3 should 
end in a triple point, where three liquid phases coexist. 
Note, that all three phase transitions occur above the 
glass transition temperature T g sa 235 K, estimated at 
zero pressure. This means, that the first transition is 
definitely a transition between two liquid phases. If T g 
does not strongly increase with pressure, the two other 
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FIG. 12: Isotherms of ST2 water from simulations in the 
restricted ensemble with 513 (open circles) and 621 (stars) 
water molecules. The density of the liquid phase at the liquid- 
vapor coexistence curve (P = bar) is shown by diamonds. 
The densities of the coexisting phases of the first and second 
liquid-liquid phase transitions of ST2 water at 260 K, obtained 
by GEMC simulations, are shown by dashed and dot-dashed 
lines, respectively. 



transitions also should be essentially liquid-liquid. 



2. ST2RF model 

Several isotherms of the ST2RF water model, obtained 
by MC simulations in the restricted NPT ensemble at T 
— 200 to 290 K, are shown in Figs. 14 and 15 together with 
the results of other simulation studies. There is good 
agreement of our data for the supercritical isotherms with 
MD simulations in the NVT ensemble (Fig. 14, T = 290 
and 275 K). In this water model the first liquid- liquid 
phase transition appears at T = 260 K and P w +1.3 
kbar in the density interval from ~0.94 to ^1.01 g cm~ 3 . 
The shape of the isotherm indicates a critical point at 
positive pressures slightly above 260 K. The shift of the 
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critical point from negative pressure for the ST2 model 
to positive pressure for the ST2RF model and the subse- 
quent disappearance of the triple point corroborates with 
the changes of the water phase diagram expected due to 
strengthening of the hydrogen bonding^ At T = 235 
K the two-phase region becomes wider (from —0.94 to 
—1.10 g cm -3 ) and the coexistence pressure increases to 
P w +2.0 kbar. Our isotherm at T = 235 K agrees well 
with the results of other simulation studies (Fig. 15). The 
systematic shift of the data points, obtained from a sed- 
imentation profile, 9 toward higher pressures is caused by 
a slightly lower temperature (230 K) and the absence of 
the LRLJ corrections in that study. 

The first liquid-liquid phase transition of the ST2RF 
model is also clearly seen in the T — 200 K isotherm, i.e. 
deeply below the glass transition temperature T g rj 255 
K, estimated at zero pressure. Probably due to the fact 
that the 200 K isotherm refers to a glassy state the dif- 
ferent branches of the isotherm start to overlap in some 
density range (Fig. 14, lower panel). This substantially 
complicates estimating the coexisting densities and the 
pressure at the transition. However, this transition re- 
mains roughly in the same pressure and density interval, 
as at T = 235 K. 

At T — 200 K the ST2RF model shows a pronounced 
second liquid-liquid (or better amorphous-amorphous) 
transition at P f=a +3 kbar with a two-phase region 
roughly from 1.1 to 1.2 g cm~ 3 . The location of the 
observed coexistence regions with respect to the liq- 
uid branch of the liquid-vapor coexistence curve of the 
ST2RF model is shown schematically in Fig. 13 (lower 
panel). Note, that the density intervals of the second 
transition of the ST2RF model and the third transition 
of the ST2 model are very close. The coexistence region 
of the first transition of the ST2RF model covers the two 
coexistence regions of the first and second transitions of 
the ST2 model. 



3. TIP4P model 

The T — 200 K isotherm of the TIP4P water model 
shows a liquid-liquid phase transition in the density in- 
terval from -1.05 to -1.08 g cm' 3 (Fig. 16). The critical 
point of this transition should be located slightly above 
200 K, at P « +1 kbar and p w 1.06 g cm" 3 . Con- 
tinuous isotherms of TIP4P water were obtained in MD 
simulations using the NVT— and NPT— ensembles (see 
stars and triangles, respectively, in Fig. 16). The shift of 
the MD data to lower densities may be attributed to the 
use of LRCI corrections in Refs0 and|2(J The compress- 
ibility maximum, obtained in the MD simulations^*^ at 
p as 1.03 g cm -3 is rather close to the critical density 
at 1.06 g cm~ 3 in our simulations. The absence of the 
liquid-liquid transition at T = 200 K in the MD simu- 
lation studies&Sii should be attributed to the limitations 
of the applied simulation methods, which underestimate 
the critical temperature (see Introduction). 
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FIG. 13: Schematic representation of the location of the 
liquid-liquid phase transitions (shadowed areas) with respect 
to the liquid branch of the liquid-vapor coexistence curve 
(solid circles) for the ST2 and ST2RF water models. Glass 
temperature at P — is indicated by dashed line. 



When decreasing the temperature to 175 K, two liquid- 
liquid phase transitions are seen. The first one is located 
at P « -0.8 kbar with densities of the coexisting phases 
of 0.96 and 1.02 g cm~ 3 . This transition is located in the 
region, which is metastable with respect to evaporation 
(see diamonds in Fig. 16). Further decrease of the tem- 
perature to T = 150 K, which is essentially below the 
glass transition temperature T g « 180 K, estimated at 
zero pressure, the different branches of the isotherm over- 
lap in a wide range of pressure and density (Fig. 16, lower 
panel). At T = 150 K the first transition apparently re- 
mains at negative pressure and its parameters (center of 
the overlap region) seem to be close to those at T = 
175 K. At T = 175 K the second liquid-liquid transi- 
tion of the TIP4P water model occurs at P « +1.8 kbar 
with densities of the coexisting phases of 1.08 and 1.14 g 
cm~ 3 (Fig. 16, middle panel). When decreasing the tem- 
perature to 150 K, this two-phase region becomes wider 
and extends up to 1.2 g cm~ 3 (Fig. 16, lower panel). 

The liquid branch of the liquid- vapor coexistence curve 
of the TIP4P model varies smoothly with temperature 
and does not show any triple point (Figs. 3, 4). Therefore, 
the first liquid-liquid transition, which occurs at negative 
pressure at T = 150 and 175 K, should end at a critical 
point between 175 and 200 K also at negative pressure. 
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FIG. 14: Isotherms of ST2RF water from MC simulations 
in the restricted NPT ensemble (open circles) and densities 
of the liquid phase at the liquid- vapor coexistence curve (di- 
amonds at P — bar). Data points from MD simulations in 
the NVT ensemble are shown by stars;- open squares- and 
solid trianglesi— 



Isotherms of TIP5P water from MC simulations in the 
restricted NPT ensemble (open circles) and densities of 
the liquid phase at the liquid- vapor coexistence curve (di- 
amonds at P = bar) . Data points from MD simulations 
in the NVT ensemble^ are shown by stars (for T = 225 
K these points were obtained by interpolation of the data 
for T = 220 and 230 K). 

The location of the two liquid-liquid transitions with 
respect to the liquid branch of the liquid-vapor coexis- 
tence curve of the TIP4P model is shown schematically 
in Fig. 17. The results of MD simulations of supercooled 
TIP4P water 6 i 23 i 26 support this picture. Indeed, the drop 
of the water density from about 1.12 to 1.05 g cm~ 3 be- 
tween 193 and 173 K along the isobar P — 2 kbar in 
Ref. 26,i, can be understood as the crossing of the second 
liquid-liquid transition of TIP4P water (Fig. 17). The 
sharp change of the water density along the isobar P « 



02 6 . probably reflects the proximity of the first liquid- 
liquid phase transition, which may be located at slightly 
positive pressures in that implementation of the TIP4P 
model (see squares in Fig. 4). 



4. TIP5P model 

Isotherms of the TIP5P water model at supercooled 
temperatures are shown in Fig. 18. The inflection point 
of the T — 250 K isotherm of the TIP5P model is lo- 
cated approximately at the same pressure and density, 
as the inflection point of the T = 200 K isotherm of 
the TIP4P model (Fig. 16). Note, that a similar shape 
of this isotherm was obtained from MD simulations in 
the NVT ensemble £ At T = 225 K a liquid-liquid phase 
transition appears at P ks +1.3 kbar with densities of 
the coexisting phases ~1.03 and ~1.08 g cm~ 3 . This 
transition was not observed at this temperature in MD 
simulations, 8 ' probably due to the unavoidable distortion 
of the subcritical isotherms in the NVT ensemble. With 
decreasing temperature the coexistence pressure of this 
transition increases noticeably and the densities of both 
coexisting phases shift to higher values. 

A lower-density liquid-liquid transition appears at T 
= 175 K with densities of the coexisting phases ^0.99 
and ~1.03 g cm~ 3 . Obviously, this transition appears 
at slightly positive pressure (see the location of the liq- 
uid density at the liquid-vapor coexistence curve at T = 
175 K indicated by a diamond in Fig. 18). At T = 150 K 
this transition is observed at negative pressures (about -3 
kbar) in the density interval from ^0.94 to ^0.99 g cm -3 . 
Thus, this liquid-liquid transition should cross the liquid- 
vapor coexistence curve at a triple point between 150 and 
175 K. Indeed, a sharp change of the liquid density at the 
liquid- vapor coexistence curve is observed at this temper- 
ature interval (see Fig. 4). The location of the estimated 
liquid-liquid two-phase regions with the respect to the 
liquid branch of the liquid- vapor coexistence curve of the 
TIP5P model is shown schematically in Fig. 17 (middle 
panel). Note some discrepancy between this phase dia- 
gram and the liquid- liquid phase transition with a critical 
point at T = 217 K, P = 3.4 kbar and p = 1.13 g cm" 3 , 
estimated for the TIP5P model in Ref. E 



5. SPCE model 

The appearance of a liquid-liquid (or better 
amorphous-amorphous) phase transition of the SPCE 
model is seen from the isotherm T = 200 K (Fig. 19). 
The critical point of this transition we estimate at 
temperature slightly above 200 K, P s=s +0.9 kbar and 
p ~ 1.06 g cm~ 3 . With decreasing temperature the 
coexistence pressure of this transition remains positive, 
whereas the density range of the two-phase coexistence 
increases from 1.07 - 1.13 g cm -3 at T — 175 K to 1.09 - 
1.16 g cm -3 at T = 150 K. A lower-density amorphous- 
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FIG. 15: Isotherm of ST2RF water at T = 235 K from 
MC simulations in the restricted NPT ensemble (open cir- 
cles) and from MD simulations in the NVT ensemble (solid 
circles-, squares-). The isotherm, obtained from a sedimen- 
tation profile at T = 230 K- is shown by stars. 

amorphous phase transition appears at slightly negative 
pressure, the coexistence region we estimate as 0.99 
- 1.04 g cm" 3 at T = 175 K and 0.97 - 1.03 g cm" 3 
at T = 150 K. The location of these two liquid-liquid 
transitions with respect to the liquid branch of the 
liquid-vapor coexistence curve is shown schematically 
in Fig. 17. Studies of the liquid-liquid transition of the 
SPCE model, based on an analytic expansion for the 
liquid free energySi do not agree with our phase diagram. 
First, the existence of a single liquid- liquid transition 
was imposed in Refl6lL The critical point of this single 
liquid- liquid transition, estimated at T « 130 K P w 
+2.9 kbar and p w 1.10 g cm~ 3 , is essentially shifted 
to lower temperature and higher pressure in comparison 
with the critical points of the two liquid-liquid phase 
transitions, shown in Fig. 17. This shift i s p robably 
caused by the use of LRCI corrections in RefroU 



IV. DISCUSSION 

Phase transitions between two liquid or amorphous 
phases with different densities were observed for sev- 
eral water models by simulations in the restricted NPT 
ensemble. The reliability of the applied method was 
tested in various ways. Firstly by comparison with direct 
GEMC simulations of the liquid-vapor and liquid-liquid 
coexistence of ST2 water. In particular, the P = iso- 
bar, simulated in the restricted NPT ensemble from the 
supercooled region up to T = 450 K well agrees with 
the liquid branch of the liquid-vapor coexistence curve 
obtained by GEMC simulations (see Fig. 2 in Ref Hoj) . 
Additionally, we have tested the sensitivity of the sim- 
ulation results in the restricted NPT ensemble to the 
number of molecules in the subcells. An increase of the 
system from 513 molecules (27 subcells, containing 19 
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FIG. 16: Isotherms of TIP4P water from MC simulations in 
the restricted NPT ensemble (open circles) and densities of 
the liquid phase at the liquid-vapor coexistence curve (dia- 
monds at P = bar). Data points from MD simulations in 
the NVT ensemble" are shown by stars and from MD sim- 
ulations in the NPT ensemble at T = 193 l& are shown 
by triangles. Experimental data points, corresponding to the 
transition between low-density and high-density amorphous 
ice at T = 110 Kp— are shown by crosses. Arrows indicate 
compression and decompression. 

molecules each) to 621 molecules (27 subcells, contain- 
ing 23 molecules each) does not affect the results notably 
(see Fig. 12, T — 260 K). This increase of the number of 
molecules in the subcells is accompanied by an increase 
of the density fluctuations in the metastable states and 
a further increase of the number of molecules per sub- 
cell can lead to phase separation in the single subcells. 
The ability of the simulations in the restricted NPT en- 
semble to locate correctly the phase transitions can be 
illustrated by using the triple (liquid-liquid-vapor) point 
of ST2 water as an example. The densities of the coex- 
isting liquid phases, obtained at about zero pressure in 
the restricted NPT ensemble are close to those obtained 
from the simulations of the liquid-vapor coexistence in 
the Gibbs ensemble and from the direct equilibration of 
two liquid phases in the Gibbs ensemble (see Fig. 12 and 
also Fig. 4 in ReflrOI). 

There is a general agreement of our simulation results 
with available simulation studies of the same water mod- 
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FIG. 17: Schematic representation of the location of the 
liquid-liquid phase transitions (shadowed areas) with respect 
to the liquid branch of the liquid-vapor coexistence curve 
(solid circles) for TIP4P, TIP5P and SPCE water models. 
Glass temperature at P ss bar is indicated by dashed line. 



els. Discrepancies are connected mainly with the differ- 
ences in implementation of the water model. In particu- 
lar, the use of long-range corrections for the intermolec- 
ular interaction in water models, which were parameter- 
ized without any long-range corrections, essentially af- 
fects the liquid density. The use of LRCI corrections 
causes a decrease of the liquid density especially at low 
temperatures. This effect is most dramatic for the ST2 
water model and therefore we consider in fact two water 
models: ST2 and ST2RF. The use of LRLJ corrections 
results in the opposite effect, i.e. in an increase of the liq- 
uid density. Again, this effect enhances with decreasing 
temperature. 

Another discrepancy between different simulation 
studies of supercooled water originates from the method, 
used for the detection (location) of the liquid-liquid phase 
transitions. Due to the intrinsic limitations of the sim- 
ulations in the canonical NVT ensemble (see Introduc- 
tion), this method notably underestimates the critical 
temperature of the phase transition. For example, in 
some temperature interval below the critical tempera- 
ture, where the simulations in the restricted ensemble 
evidence a phase transition, the isotherms, obtained in 
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FIG. 18: Isotherms of TIP5P water from MC simulations in 
the restricted NPT ensemble (open circles) and densities of 
the liquid phase at the liquid-vapor coexistence curve (dia- 
monds at P = bar). Data points from MD simulations in 
the NVT ensemble- are shown by stars (for T = 225 K these 
points were obtained by interpolation of the data for T = 220 
and 230 K). 



the NVT ensemble, have a supercritical appearance (see 
Fig.16, T = 200 K; Fig.18, T = 225 K; Fig.19, T = 
200 K). Note also, that the applicability of simulations 
in the NVT ensemble for the study of phase transitions is 
doubtful in the case of multiple phase transitions, which 
occur in a narrow pressure range (see, for example Fig. 12, 
T = 235 and 260 K), due to unavoidable and unpre- 
dictable distortions of the subcritical isotherms. 

The phase diagrams of the liquid- vapor and the liquid- 
liquid phase transitions are highly sensitive to the wa- 
ter model. All studied water models (except the ST2RF 
model) show comparable liquid densities at ambient con- 
ditions, where all these models were parameterized. Both 
upon heating and cooling, these models show increas- 
ing deviations from the saturated liquid density. For in- 
stance, the liquid-vapor critical temperature T c varies 
from ~ 540 K to ~ 640 K for the different water models 



(see Table I). The influence of the water model is even 
more obvious, when the distance between T c and the 
temperature T max of the liquid density maximum is con- 
sidered (Table I). This distance is largest and comparable 
with the experimental value for the three-site SPCE wa- 
ter model, whereas its is much smaller for five-sites water 
models (ST2, ST2RF and TIP5P). Below T max and at 
supercooled temperatures the liquid density at the liquid- 
vapor coexistence curve and the densities of the coexist- 
ing liquid phases of the liquid- liquid phase transitions are 
highly sensitive to the water model and its implementa- 
tion. 

The multiplicity of liquid-liquid phase transitions, 
which we obtained for five non-polarizable models and 
Jedlovszky and Vallaur: 12 obtained for a polarizable wa- 
ter model, seems to be more comprehensible, than the 
existence of only one liquid-liquid transition. Obviously, 
these transitions are related to the variety of amorphous 
phases of the one-component isotropic fluid, which dif- 
fer in their local order. It is natural to assume, that 
each crystalline phase (or at least the crystalline phases, 
which can directly transform to the liquid phase) should 
result in a corresponding amorphous (liquid) phase with a 
short-range order, which is reminiscent of the crystalline 
phase's So, a multiplicity of liquid-liquid phase tran- 
sitions may be expected for water, showing at least 12 
crystalline phases, 5 of which are in direct equilibrium 
with the liquid'. This corroborates with the existence 
of at least 3 amorphous phases of supercooled water: 
low-density amorphous ice (LDA), high-density amor- 
phous ice (HDA)i and very-high-density amorphous ice 
(VHDA)' Moreover, the existence of other amorphous 
phases of supercooled water can not be excluded. 68 '- 9 

The existence of more than one liquid-liquid phase 
transition in a one-component system was also con- 
cluded from simulation studies of other model fluids'^ 
An isotropic soft-core intermolecular potential with two 
steps causes the appearance of one liquid-liquid phase 
transitions^ Two liquid-liquid phase transitions were 
observed, when a third step was included in the soft-core 
potential, 7 - Two liquid-liquid phase transitions were also 
observed for a waterlike model fluid confined between hy- 
drophobic surfaces' 

We have found three liquid-liquid phase transitions of 
the ST2 water model. Hence, there are 4 distinct phases 
of supercooled water, which we have numbered as phase 
I to IV, starting from the phase with the lowest den- 
sity. The density range, where the liquid-liquid tran- 
sitions occur, is found to be the widest for ST2 and 
ST2RF water (this coincides with the wide hysteresis 
loop obtained for the ST2RF model by quick isother- 
mal compression/decompression 22 ). The density interval 
of the two lower-density transitions in the ST2 model 
practically coincides with the density range of the first 
transition in the ST2RF model. Thus, 3 supercooled 
phases (phases I, III and IV) are found for the ST2RF 
model. Two liquid-liquid phase transitions and 3 super- 
cooled phases (phases I and, presumably, phases II and 
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FIG. 19: Isotherms of SPCE water from MC simulations in 
the restricted NPT ensemble (open circles) and densities of 
the liquid phase at the liquid-vapor coexistence curve (dia- 
monds at P = bar). Data points at T = 200 K from MD 
simulations in the NVT ensemble— are shown by stars. 



Ill) we have found for TIP4P, TIP5P and SPCE water. 
The density range of the liquid-liquid transitions in each 
of these models is much narrower than in the ST2 and 
ST2RF models. The latter observation corroborates the 
narrow hysteresis loop obtained for the TIP4P model by 
quick isothermal compression/decompression. 22 

The density of LDA (the lowest-density phase of amor- 
phous water observed experimentally) is about 0.94 g 
cm~ 3 at zero pressure and T = 11 Ki At the same condi- 
tions, the density of HDA was reported from 1.14-1.15 g 
cm -3 in Refl67lto 1.17-1.19 g cm -3 in 

Ref[l|) , 

whereas the 

density of VHDA is about 1.25-1.26 g cm~ 3 , 67 Comparing 
with the densities of our isotherms at the same pressure 
we can attribute LDA to the lowest density phase (phase 
I) of the first liquid-liquid phase transition of all studied 
water models. HDA likely corresponds to phase III in 
all studied models, whereas VHDA could be attributed 
to phase IV, which is observed in ST2 and ST2RF water 
only. 

Phase II, which is observed for all studied water mod- 
els, except the ST2RF model, usually corresponds to liq- 
uid water in equilibrium with the vapor (see the dia- 
monds in Figs. 12, 16, 18 and 19): in the whole studied 
temperature range of the TIP4P and SPCE water mod- 
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els, above the triple point for the ST2 model and be- 
low the triple point for the TIP5P model. The den- 
sity of phase II in the supercooled region of the ST2, 
SPCE and TIP5P models is about 1.00-1.04 g cm~ 3 and 
varies slightly with the temperature. For the TIP4P wa- 
ter model phase II achieves the density 1.14 g cm -3 with 
decreasing temperature at zero pressure. The search of 
such a phase in real supercooled water seems to be a 
challenge for the experimental studies. There has been a 
single experimental observation of amorphous water (hy- 
perquenched water) with the density 1.04 g cm -3 , re- 
ported in Reflrisl 

Experimentally, transitions between various amor- 
phous ices can be studied directly at temperatures below 
T w 150 K, i.e. below the temperature of the sponta- 
neous crystallization of the amorphous water. At these 
temperatures the phase transition between HDA and 
LDA occurs at a pressure of P s» 2 kbar^SSiZi which 
slightly increases upon coolingi 7 * Measurements of the 
low-temperature metastable melting curves of crystalline 
icesS provide an extension of the HDA/LDA transition 
line up to about 230 K and down to pressures of about 
+0.5 kbar. The densities of the coexisting phases at 
the HDA/LDA phase transition is about 0.94 and 1.20 
g cm~ 3 at T = 130 and they do not change essen- 
tially upon cooling to T = 110 K. 62,63 Experimental data 
points for the isotherm T = 110 K^iS 3 . are compared 
with the simulated T — 150 K isotherm of TIP4P wa- 
ter in Fig. 16. The experimental hysteresis loop approxi- 
mately covers the density interval, which corresponds to 
the two liquid-liquid (amorphous-amorphous) transitions 
of the TIP4P, TIP5P and SPCE water models. This con- 
firms our assignment of phase I to LDA and phase III to 
HDA. A comparison of the experimental and simulated 
data, shown in Fig. 16 (lower panel) may also be used as 
an indication for the existence of phase II in real water. 
The location of the phase transitions of the TIP4P model 
in the pressure-density plane is closer to the experiment, 
than of TIP5P and SPCE models. This conclusion is in 
line with the results of computer simulations^**! 7 , which 
show the ability of the TIP4P model to describe cor- 
rectly the phase diagram of crystalline ices, unlike the 
SPCE and TIP5P models. Note finally, that the ST2 
water model likely underestimates the densities of HDA 
and VHDA due to the overestimation of the tetrahedral 
water structure. 

Other common (universal) features of the phase be- 
havior can be noticed. In particular, the second transi- 
tion in all studied models (except the specific case of the 
ST2RF model) seems to be rather universal. It is always 
located at positive pressures and its coexistence pressure 
increases with decreasing temperature. The proximity 
of the critical point of this transition to the liquid den- 
sity maximum and to the consecutive density minimum 
suggests to attribute the existence of these features to 
the influence of the second liquid-liquid phase transition. 
Indeed, above T max the influence of the liquid-vapor crit- 
ical point on the density of the liquid branch is dominant 



and the temperature derivative of its density is negative, 
i.e. it shows normal behavior. Below T rnax the deriva- 
tive (dp/dT)p is positive and the liquid branch can be 
considered as a supercritical isobar (P « 0) of the sec- 
ond liquid-liquid phase transition. With further decreas- 
ing temperature, the zero pressure isobar turns to the 
"normal" behavior with (dp/dT)p < after the density 
minimum. This "normal" behavior of the isobar (P w 
0) in deeply supercooled region can be attributed to the 
tilt of the two-phase region of the second liquid-liquid 
phase transition toward higher densities with decreasing 
temperature (Fig. 17) rather than to the influence of the 
distant liquid-vapor critical point. So, the liquid density 
minimum and its consecutive maximum appear to be the 
result of a crossover between these three regimes. Note, 
that the critical temperature of the second liquid-liquid 
phase transition in all studied models appears above (or 
close to) the glass transition temperature. 

The first liquid-liquid phase transition is much more 
sensitive to the choice of the water model and its imple- 
mentation. This transition is located completely at posi- 
tive pressures only for the ST2RF water model, which 
does not provide a realistic water density at ambient 
pressures. It is located completely at negative pressures 
(TIP4P and SPCE models), or it meets the liquid- vapor 
coexistence curve in a triple point (ST2 and TIP5P mod- 
els). The critical point of this transition is located at 
slightly negative pressures for the ST2, TIP4P and SPCE 
models and at slightly positive pressure in the case of the 
TIP5P model. Only for the ST2 and TIP4P models the 
critical point of the first liquid- liquid phase transition is 
above the estimated zero pressure glass transition tem- 
perature. For the TIP5P and SPCE water models this 
transition has be considered as a transition between two 
glassy states. 

The line of the first liquid-liquid phase transition in 
the P-T plane is characterized by an essentially nega- 
tive slope for the ST2 model and an essentially positive 
slope for the TIP5P model, whereas it could not be de- 
termined for the TIP4P and SPCE models. The presence 
of a triple point, where the vapor coexists with two liq- 
uid phases, means, that the liquid branch of the liquid- 
vapor coexistence curve crosses the liquid-liquid spinodal 
at some temperature below the triple point temperature. 
Such crossing can explain the apparent singular behav- 
ior of some thermodynamic properties, which in real liq- 
uid water is observed at w 228 K, 78 i.e. 49° below the 
temperature of the density maximum. This scenario is 
very close to that observed for the ST2 model, where the 
triple point is located about 35° below T max and the 
liquid-liquid spinodal crosses the liquid- vapor coexistence 
curve at about 40° below T max m& The triple point for 
the TIP5P model is found about 100-125° below T max 
and therefore the liquid-liquid spinodal should cross the 
liquid- vapor coexistence curve far away from T max . 

The variation of the heat capacity of liquid water with 
decreasing temperature seems to be highly sensitive to 
the location of the first liquid-liquid phase transition. In 
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the TIP4P and SPCE water models, which do not show a 
triple point, noticeable temperature changes of C p upon 
cooling can not be detected. A pronounced increase of 
C p upon cooling is observed for the TIP5P, ST2 and 
ST2RF water models. The first two models show a triple 
point of the first liquid- liquid and the liquid- vapor phase 
transition. 

The structure of liquid water in various thermody- 
namic states can be studied by the analysis of cluster- 
ing and percolation^ In particular, the apparent sin- 
gular behavior of some properties of liquid water in the 
supercooled region was attributed to the percolation of 
"four-bonded" water molecules J2, This idea can be in- 
corporated naturally in the observed phase behavior of 
supercooled water. Indeed, the percolation line of physi- 
cal clusters of one of the constituting components should 
coincide with the spinodal of the phase transition, cross 
the coexistence curve in the critical point and extend 
into the one phase regioni^iiS Thus, crossing the spin- 
odal of a phase transition is equivalent to crossing a per- 
colation line. The analysis of water clustering allowed 
the localization of the percolation line in aqueous solu- 
tions with respect to the liquid-liquid phase transition. 83 
We have applied this method also to identify the molec- 
ular species with different local order, which may be re- 
sponsible for the liquid-liquid phase transitions in super- 
cooled ST2 waterii2i2i The first, lowest-density liquid- 
liquid transition is caused by the formation of an infinite 
hydrogen-bonded network of perfectly coordinated tetra- 
hedral water molecules. With decreasing temperature 
along the liquid-vapor coexistence curve, the spinodal of 
the first liquid-liquid phase transition is crossed, which is 
the percolation line of the perfectlycoordinated tetrahe- 
dral water molecules (Fig. 4 in Ref 10). 

Based on our simulations of the phase diagrams of vari- 
ous water models in the supercooled region, we can spec- 
ulate about the possible location of liquid-liquid phase 
transitions in real water. We may expect, that there are 
at least two liquid-liquid phase transitions in real water, 
which effect the properties of liquid water at ambient 
conditions. One transition should be located at posi- 
tive pressures with its critical point close to the density 
maximum of liquid water. This transition seems to be re- 
sponsible for the anomalous behavior of the liquid water 
density. Another transition may meet the liquid-vapor 
transition and the spinodal of this transition could be 
the origin of the experimentally observed thermodynamic 
singularities of liquid water at 228 K. 78 

Finally, we would like to mention some perspectives 
for further studies of the phase behavior of fluids. The 
existence of various liquid phases of a one-component 
isotropic fluid is related presumably to different kinds of 
short and intermediate range order in the liquids. Indeed, 
the liquid phase with the lowest density consist mainly 
of perfectly tetrahedrally coordinated water molecules, 
whereas the liquid water phase with the highest den- 
sity consists mainly of water molecules with close packed 
Lennard-Jones-like local order^fl The main structural 



features of the other water phase(s) remain unclear and 
deserve further studies. The structure of these phases 
has to be analyzed in more detail, taking into account 
the various kinds of local ordering in crystalline ices. 

In accord with the liquid-vapor phase transition, the 
fluid density is usually considered as the order parameter 
of the liquid-liquid phase transition of a one-component 
fluid. On the other hand, the concentration is the or- 
der parameter of liquid-liquid phase transitions in mix- 
tures. In the case of a one-component fluid it seems to 
be reasonable to consider the concentration of molecules, 
showing some definite kind of local order, as an order 
parameter of the liquid-liquid phase transition. Such 
an approach is reminiscent of the well-known two- and 
multi-states water models (see, for example Ref 85), orig- 
inating from the ideas of HWhiting and W.C.Roentgen 
of the 19th century*^ 7 . Note, that the mutual trans- 
formations between different kinds of local ordering in a 
one-component fluid violates the conservation of the con- 
centrations of species. Thus, finding a more appropriate 
order parameter of the liquid-liquid phase transition in a 
one-component fluid is an actual problem of the physics 
of liquids. Moreover, the universality class of this transi- 
tion is not necessarily the universality class of the Ising 
model, as in the case of binary mixtures. For example, 
the possibility to change the local order continuously in- 
troduces unavoidable disorder, which could vary with the 
thermodynamic conditions. This means, that the critical 
behavior of such a system could belong to the universality 
class of random field Ising models. We can not exclude, 
that the enhancement of disorder with increasing tem- 
perature could lead to a rounding of the phase transition 
and the disappearance of a true critical point. In this 
case the two coexisting liquid phases are not infinite and 
the two-phase state appears as a domain structure. 
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